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1. Introduction 



This paper deals with the pricing of options depending on several underlying assets. While there is 
a vast amount of literature on the pricing of single-asset options, see, e.g., [9, 42] for an overview, 
the amount of literature considering the multi-asset case is rather limited. This is most likely due to 
the fact that the trade-off between flexibility and tractability is particularly delicate in a multivariate 
setting. On the one hand, the model under consideration should be flexible enough to recapture 
stylized facts observed in real option prices. When dealing with multiple underlyings, this becomes 
challenging, since not only the individual assets but also their joint behaviour has to be taken into 
account. On the other hand, one needs enough mathematical structure to calculate option prices in 
the first place and to be able to calibrate the model to market prices. Due to an increasing number of 
state variables and parameters, this is also not an easy task in a multidimensional framework. In this 
article we propose the multivariate OU-type stochastic volatility model of Pigorsch and Stelzer [38] 
in the generalised form introduced by Barndorff-Nielsen and Stelzer [4], which seems to present a 
reasonable compromise between these competing requirements. 

The log-price processes Y = (Y l , . . . ,Y d ) of d financial assets are modelled as 



where ii E M. d , A is a real d xd matrix, and j3, p are linear operators from the real d x d matrices to 
M. d . Moreover, W is an M^-valued Wiener process and L is an independent matrix subordinator, i.e., a 
Levy process which only has positive semidefinite increments. Hence, the covariance process £ is an 
Ornstein-Uhlenbeck (henceforth OU) type process with values in the positive semidefinite matrices, 
cf. Barndorff-Nielsen and Stelzer [3]. Thus we call (1.1), (1.2) the multivariate stochastic volatility 
model of OU type. The positive semidefinite OU type process £ introduces a stochastic volatility 
and, what is difficult to achieve using several univariate models, a stochastic correlation between the 
assets. Moreover, £ is mean reverting and increases only by jumps. The jumps represent the arrival 
of new information that results in positive shocks in the volatility and positive or negative shocks 
in the correlation of some assets. Due to the leverage term p{dL t ) they are correlated with price 
jumps. The present model is a multivariate generalisation of the non-Gaussian OU type stochastic 
volatility model introduced by Barndorff-Nielsen and Shephard [2] (henceforth BNS model). For one 
underlying, these models are found to be both flexible and tractable in Nicolato and Venardos [37]. 
The key reason is that the characteristic function of the return process can often be computed in closed 
form, which allows European options to be be priced efficiently using the Fourier methods introduced 
by Carr and Madan [8] as well as Raible [39]. In the present study, we show that a similar approach 
is also applicable in the multivariate case. Recently, Benth and Vos [5] discussed a somewhat similar 
model in the context of energy markets. However, they do not establish conditions for the applicability 
of Fourier pricing and, more importantly, do not calibrate their model to market prices. 

Alternatively, the covariance process £ can also be modelled by other processes taking values in the 
positive semidefinite matrices. In particular, several authors have advocated to use a diffusion model 
based on the Wishart process, cf., e.g., Da Fonseca, Grasselli, and Tebaldi [13], Gourieroux [20], 
Gourieroux and Sufana [21], and Da Fonseca and Grasselli [11]. This leads to a multivariate gener- 
alisation of the model of Heston [24]. However, there is empirical evidence suggesting that volatility 
jumps (together with the stock price), cf. Jacod and Todorov [31], which cannot be recaptured by a 
diffusion model. Moreover, the treatment of square-root processes on the cone of positive semidefinite 



dY t = 
dLt = 



(il + P %))dt + Zj dW t + p (dLt) 
(AI lt + Z t A T )dt + dL t , 



(1.1) 
(1.2) 



2 



Option Pricing in Multivariate Stochastic Volatility Models of OU Type 



matrices is mathematically quite involved, see Cuchiero, Filipovic, Mayerhofer, and Teichmann[10]. 
For example, whereas Da Fonseca and Grasselli [11] have very recently succeeded in calibrating their 
model to market prices, the resulting parameters do no satisfy the drift condition for the existence 
of the underlying square-root diffusion, suggesting that a more sophisticated optimization routine is 
necessary. 

Another possible approach is to consider multivariate models based on a concatenation of univariate 
building blocks. This approach is taken, e.g., by Luciano and Schoutens [34] using Levy processes, 
by Dimitroff, Lorenz, and Szimayer [14], who consider a multivariate Heston model, and by Hubalek 
and Nicolato [27], who put forward a multifactor BNS model. However, all these models either 
have a somewhat limited capability to catch complex dependence structures (compare Section 4.2) or 
lead to tricky (factor) identification issues. Apart from models where all parameters are determined by 
single-asset options, we are not aware of successful calibrations of such models. The paper of Ma [35] 
proposes a two-dimensional Black-Scholes model where the correlation between the two Brownian 
motions is stochastic and given by a diffusion process with values in an interval contained in [—1,1]. 
However, pricing can only be done via Monte-Carlo simulation in this model. In addition, an extension 
to higher dimensions is not obvious, since the necessary positive semidefiniteness of the correlation 
matrix of the Brownian noise imposes additional constraints, which are hard to incorporate. 

The remainder of this paper is organised as follows. Sections 2.1 and 2.2 introduce the multivariate 
stochastic volatility model of OU type. Afterwards, we derive the joint characteristic function of 
(Y t ,L t ). We then show in Section 2.4 that a simple moment condition on L implies analyticity and 
absolute integrability of the moment generating function of Y t in some open complex strip around 
zero. Equivalent martingale measures are discussed in Section 2.5, where we also present a subclass 
that preserves the structure of our model. In Section 3, we recall how to use Fourier methods to 
compute prices of multi-asset options efficiently. Subsequently, we propose the OU-Wishart model, 
where L is a compound Poisson process with Wishart distributed jumps. It turns out that the OU- 
Wishart model has margins which are in distribution equivalent to a T-OU BNS model, one of the 
tractable specifications commonly used in the univariate case. Moreover, the characteristic function 
can be computed in closed form, which makes option pricing and calibration particularly feasible. 
In an illustrative example we calibrate a bivariate OU-Wishart model to market prices, and compare 
its performance to the multivariate Variance Gamma model of [34] and a multivariate extension with 
stochastic volatility. As a final application, we show in Section 5 that covariance swaps can also be 
priced in closed form. The appendix contains a result on multidimensional analytic functions which 
is needed to establish the regularity of the moment generating function in Section 2.4. 

Notation 

Mj „(E) (resp. M^„(C)) represent the d x n matrices with real (resp. complex) entries. We abbreviate 
Mrf(-) = M c i.d{-). §d denotes the subspace of M^(R) of all symmetric matrices. We write §^ for the 
cone of all positive semidefinite matrices, and §j~ + for the open cone of all positive definite matrices. 
The identity matrix in Mj(M.) is denoted by I c /. a (A) denotes the set of all eigenvalues of A £ Mj(C). 
We write Re(z) and Im(z) for the real or imaginary part of z £ C d or z £ M^(C), which has to be 
understood componentwise. The components of a vector or matrix are denoted by subscripts, however 
for stochastic processes we use superscripts to avoid double indices. 



This study generalizes the theory of affine processes from the positive univariate factors treated in [16, 15] to factor 
processes taking values in the cone of symmetric positive semidefinite matrices. In particular, to ensure the existence of 
square-root processes, a quite intricate drift condition turns out to be necessary. 
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On W l , we typically use the Euclidean scalar product, (x, y) Rd := x T y, and on Mj(M.) or the 
scalar products given by {A,B) Md ^ := tr(A T S) or (A,B) Sd := tr(AB), respectively. However, due to 
the equivalence of all norms on finite dimensional vector spaces, most results hold independently of 
the norm. We also write (x,y) = x T y for x, y £ C d , although this is only a bilinear form but not a scalar 
product on C rf . 

We denote by vec : M ( j(R) — > M. d2 the bijective linear operator that stacks the columns of a matrix 
below one other. With the above norms, vec is a Hilbert space isometry. Likewise, for a symmetric 
matrix S £ §^ we denote by vech(5) the vector consisting of the columns of the upper-diagonal part 
including the diagonal. 

Furthermore, we employ an intuitive notation concerning integration with respect to matrix- valued 
processes. For an M„ V ,(IR)- valued Levy process L, andM^ m (]R) resp. M n . p (R)- valued processes X,Y 
integrable with respect to L, the term JqX s dL s Y s is to be understood as the d x p (random) matrix with 
(/, 7 >th entry £f =1 £"=i &X* dL>?Y l s j . 

2. The multivariate stochastic volatility model of OU type 

For the remainder of the paper, fix a filtered probability space {£1,^ ', {^t)te[o,T]^P) i n th e sense of 
[30, Definition 1.1.3], where = {^0} is trivial and T > is a a fixed terminal time. 

2.1. Positive semidefinite processes of OU type 

To formulate our model, we need to introduce the concept of matrix subordinators as studied in [1]. 

Definition 2.1. An ^-valued Levy Process L = (L f is called matrix subordinated-, ifL t — L s £ 
for all t > s. 

The characteristic function of a matrix subordinator L is given by £(e /tr ( ZL ')) = exp(> L (Z)) for the 
characteristic exponent 

Wl(Z) = HylZ) + [ (e h < xz *> - 1) K L (dX), Z e M d (R), 
where Yl £ §J and Ki is a Levy measure on satisfying fO^S^S^) = as well as ||x||<i} 1 1^1 1 K t{dX) < 

oo. 

Positive semidefinite processes of OU type are a generalisation of nonnegative OU type processes 
(cf. [3]). Let L be a matrix subordinator and A £ M</(]R). The positive semidefinite OU type process 
£ = (r ? )(eR + is defined as the unique strong solution to the stochastic differential equation 

dZ t = (AL t + L t A T )dt + dL t , L Q £§+. (2.1) 

It is given by 

L t = e A %e AT ' + f dL s e^'-*) . (2.2) 

Jo 

Since £ f £ for all t £ E + , this process can be used to model the stochastic evolution of a covariance 
matrix. As in the univariate case there exists a closed form expression for the integrated volatility. 
Suppose 

0£o(A) + o(A). (2.3) 
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Then, the integrated OU type process £ + is given by 

£+ := f L s ds = A- 1 (L t -Iq-L,), (2.4) 
Jo 

where A:X4 AX +XA T . Note that condition (2.3) implies that the operator A is invertible, cf. 
[26, Theorem 4.4.5]. In the case where £ is mean reverting, i.e., A only has eigenvalues with strictly 
negative real part, condition (2.3) is trivially satisfied. 

2.2. Definition and marginal dynamics of the model 

The following model was introduced and studied in [38] from a statistical point of view in the no- 
leverage case and has also been considered in [4]. Here we discuss its applicability to option pricing. 

Let L be a matrix subordinator with characteristic exponent 1//^ and W an independent Revalued 
Wiener process. The multivariate stochastic volatility model of OU type is then given by 

dY t = (ji + Pfa))dt+I?dW t + p(dL,), Y Q eR d (2.5) 
dZ, = (AL t + L t A T )dt + dL t , LoeSj, (2.6) 

with linear operators p\p : M d (R) ->■ W , p G R d , and A G M d (R) such that £ a (A) + a (A). 

We have specified the risk premium jS and the leverage operator p in a quite general form. The 
following specification turns out to be particularly tractable. 

Definition 2.2. We call j8 and p diagonal if, for jSi , . . . , j5 c i G R and Pi, . • -,pd G M, 

/ jSiXn \ / piXn 



p(X)=\ : |, VXGM rf (R). 



In the following, we will denote for each i G {1, . . . , d} by j3 ! (X) and p ! (X) the z'-fh component 
of the vector j8(X) or p(X), respectively. The marginal dynamics of the individual assets have been 
derived in [4, Proposition 4.3]: 

Theorem 2.3. Let i G {1, . . . ,d}. Then we have 

\ Jo J teR+ 

where ^= denotes equality of all finite dimensional distributions. 

fax \ 

Let us now consider the case where A is a diagonal matrix, A = ■ . , and j8 , p are diagonal 



as well. Then, for every i G {1, ...,</}, we have 



£0? 7 = (H + ^dt + ffiidWi + PidLf, (2.7) 
c/lf = la^dt + dV}. (2.8) 

Evidently, every diagonal element L", i = 1, . . . ,d, of a matrix subordinator L is a univariate subor- 
dinator, and thus £" is a nonnegative OU type process. Consequently, the model for the z'-th asset is 
equivalent in distribution to a univariate BNS model. 
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2.3. Characteristic function 



Let (•, (v)w b e bilinear forms as introduced in the notation, where V,W may be either M. d , C d or 
M c /(-). Given a linear operator T : V — > W, the adjoint T* : W — > V is the unique linear operator such 
that (Tx,y} w = (x, T*y) v for all x G V and )i£ff. Directly by definition we obtain the following: 

Lemma 2.4. Let y G W 1 , z G Mj(R) cotg? ? G R + . Then the adjoints of the linear operators 

A:X ^AX + XA T , m(t) : X \-+ e A, Xe Al> -X, 
tf(t) : X h+ AA + j3 {A-\&{t)X))y T + p{X)y T + ^/A" 1 

on M c i(C) are given by 

A*:X^A T X + XA, £g(t)* :X ^ e AJ 'Xe At -X, 
V(tY :X^e AT, Xz T e At +p*(Xy)+^(t)*A-* (p*(Xy)+ l -XyyA . 



Note that for diagonal p it holds that p*(X) 



for all X G M c 



\ o p d x dd / 

Our main objective in this section is to compute the joint characteristic function of (Y t ,L t ). This will 
pave the way for Fourier pricing of multi-asset options later on. Note that we use the scalar product 



{(x h y l ),(x 2 ,y 2 )) ■=xJx 2 +tT(yJy 2 ) 



on W 1 xM d (R). 



Theorem 2.5 (Joint characteristic function). For every (y,z) G R d xM d (R) and t G R+, the joint 
characteristic function of(Y t ,'L t ) is given by 

E[exp(i((y,z),(Y t ,Lt)))] 
= exp { iy T (Y + pt) + itr(E e ATf »? A ') 

+/tr (lo (V Tf A-* U*{y) + l -yyA e At - A~* U*(y) + ^yy T 

+ jf Yfc (e^ V + p* (y) + e^^A-* U* (y) + lyyA e As - A"* (V (y) + ^yy T ) ) ds 

where A - * := (A*) -1 denotes the inverse of the adjoint of A : X i— > AX +XA T , f/jaf zj, inverse of 
A* :X^A T X+XA. 

Note that for z = we obtain the characteristic function of Y t . 

Proof. Since £ is adapted to the filtration generated by L, and by the independence of L and W, 



E[exp(((y,z),(Y t ,L t )))] = e iyT{Yo+ ^E 
= e iy J (Yo+nt) E 



exp( --y Ify 
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By (2.4) and using the fact that the trace is invariant under cyclic permutations the last term equals 

e iy^(Y +Ht) E l e itr(z T I.,+l}(A-\Z l -Zo-L,)) } ^+p(L,)y^+^A- l (Z,-L -L,)i _ 

In view of (2.2), we have 

E f - Eo -Lt = f 2%{t - s) dh s + @(t)Lo, 
Jo 

for the linear operator 38{t) from Lemma 2.4. Therefore, 

E[exp(i((y,z),(Y t ,L t ))} 

= exp (fy T (Y + lit) + /tr f z Vlo^' + J8 (A" 1 ^/)^))/ + ^A" 1 (#(*)*<>)) ) 
x E exp (/tr (V J\ A ^ dL s e A '^ + J8 (a" 1 (^JV^Ul^^ -L^jy 1 
+ P (Lt)y T + ixy T A- 1 ( jf ^ ^ ^ T (M - L,) ) )" 
= exp (iy T (Y + iit) + it r (zV'E /' + /J(A- 1 (#(f)Io))y T + p^A"^^) 



x £ 



exphtrf (£tf(t-s)dL^J I d 



with the linear operator ^{t) from Lemma 2.4, since A 1 y | '^ ( ' a) dL s e^^ ^ - L,j £ § d . An im- 
mediate multivariate generalisation of results obtained in [40, Proposition 2.4] (see also [18, Lemma 
3.1]) yields an explicit formula for the expectation above: 



exp /tr 



J^(t-s)dL^ I d Jj = e xp(j\ L (tf(s)*I d )ds\ 



By Lemma 2.4 we have 

e S' VL (V( S yi d )ds = Jo YL(e AT h T e A '+p*(y)+e AT -'A-*(p*(y)+^y T )e As -A-*(p*(y)+^yy T ))d S 

This expression is well-defined, because 

e ATs z T e As + p*(y)+e ATs A-* U* ( y ) + l - yy A <A S - U*{ y )+ l - yy A eM d (R) + i§+ 



d ■ 



for all s G [0,/]. Indeed, this follows from 

e^A-*(^ T )^-A-*(yy 
Finally, we infer from Lemma 2.4 that 



e ATu yy T e Al 'dueS+ 



(2.9) 



tr (^(A-'^Olo))/ + ^yy T A- l m)L )j = tr ( £ A"* \B*(y) + l -yy J 

which gives the desired result by noting that tr(z£ f ) = tr(z T £ ( ). 



□ 
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2.4. Regularity of the moment generating function 

In this section we provide conditions ensuring that the characteristic function of Y t admits an ana- 
lytic extension <I>y ( to some open convex neighbourhood of in C d '. Afterwards, we show absolute 
integrability. The regularity results obtained in this section will allow us to apply Fourier methods in 
Section 3 to compute option prices efficiently. 

Definition 2.6. For any t G [0, T], the moment generating function ofY t is defined as 

<S> Yt {y) :=E[fssp(y T Y t )], 

for all y S C d such that the expectation exists. 

Note that <J>y ( may not exist anywhere but on iM. , where it coincides with the characteristic function 
of Y t . The next lemma is a first step towards conditions for the existence and analyticity of the moment 
generating function <t>y, in a complex neighbourhood of zero. 

Lemma 2.7. Let Lbe a matrix subordinator with cumulant transform &l, that is 

& L (Z) = Yd-iZ) = tr( 7L Z) + / (e a ^ - 1) K L (dX), Z £ M d (C), 

and let s > 0. Then ®l is analytic on the open convex set 

S e :={ZEM d (C):\\Re(Z)\\ <£}-§+, (2.10) 

if and only if 

f e tr(RX) KL ( d x) < oo for all Re M d (R) with \\R\ \ < e. (2.11) 

7{||x||>i} 

Proof If (2.11) holds, [15, Lemma A.2] implies that Z h-> E(e t ' c{ZL ^) = e &L ^ is analytic on S e . Due 
to Assumption (2.11), dominated convergence yields that 0^ is continuous on S e . The claim now 
follows from Lemma A.l. Conversely, if 0^ is analytic on S e , then [15, Lemma A.4] implies that 
£(tr(zz4)) = e & L (z) for all zeS E . Thus, by [41, Theorem 25. 17], Condition (2.1 1) holds. □ 

The next theorem is a nontrivial (especially due to the involved heavy matrix calculus) generaliz- 
ation of [37, Theorem 2.2] to the multivariate case. It holds for all sub-multiplicative matrix norms 
on Md(M) that satisfy | |_y_y T 1 1 = | \y\ | 2 for all y £ W 1 , where we use the Euclidean norm on M. d . For 
example, this holds true for the Frobenius and the spectral norm (the operator norm associated to the 
Euclidean norm). 

Theorem 2.8 (Strip of analyticity). Suppose the matrix subordinator L satisfies 

f e u(RX) Kh ( dX ) < oo for all R G M d (M) with I \R\ I < e, (2. 12) 

J{\\X\\>i} 

for some £ > 0. Then the moment generating function <£>y ( ofY t is analytic on the open convex set 

S e :={yeC:||Re(y)||<e}, 

where 

8:= - (^fijiiA--ir | "" |+vs>0 (2 - 13) 
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+ P + 



2s 



Moreover, 



.(gW + l^lA- 1 !! " r "J (e 2 \\ A W> + l)\\A- l \\ 

® Yt (y)=exp ^y T (Yo + ^)+t r (L Q H y (t)) + ^e L (H y (s)+p*(y)) ds^j (2.14) 
for all y G Sq, where 

H y (s) := e ATs A-* (f3*(y) + lyyA - A"* U*{y) + lyyA . (2.15) 

Proof. The main part of the proof is to show that the function 

G(y) :=exp (y T (Y + p.t) +tr(Eo# y (*)) + £® L (H y (s) + p*(y)) ds 



is analytic on Sq. First we want to find a 8 such that for all u G R^ with < 8, it holds that 

\\H u (s) + p* (u)\\ < e for all s G [0,f]. Since 



<I(,2||A|| f+1) || A -l|| ||M|| 2 + ^ ||p|| + (e 2||A| k + 1) || A -l||^ 



+p*(«) 



we have to find the roots of the polynomial 

^^^(^^^^(^^^^(lIpll + ^l^^ + ^llA^lll^l^x-s. 

The positive one is given by 8 as stated in (2.13). Note that 8 > 0, because p is a cup-shaped parabola 
with p(0) = -e < 0. 

Now let y £ Sq, i.e., y = u + iv with | |u| | < 8. Using Re(yy T ) = uu T — vv T and (2.9) we get 

Re(H y (s)+p*(y)) = H u (s)+p*(u) - \ (e ATs A-*(vv T )e As - A-*(vv T )) 

= H u (s)+p*(u) - 1 fe ATr vv T e Ar dr. 
2 jo 

Because of /q e A r vv T e Ar dr G Sj, we have 

/ e tr(Re(// T (.v)+p*(v))X) 
■/{||X||>U 

e tr((//„(.v)+p*( U ))X) e -2"VV J0,; " c U 'D K L {dX) <oo 



'{||*II>1} 

w tr((^(5)+p*(»))X) e -3t r ((/o« ATr vv T e Ar <ir)x 
'{PH>1}" 

by Assumption (2.12), since | \H u (s) + p*(u)|| < £■ Thus, by Lemma 2.7 the function 

S e Gy^0 L (// y (*)+p*(y)) 

is analytic on Sq for every s G [0,?]. An application of Fubini's and Morera's theorem shows that 
integration over [0,?] preserves analyticity, cf. [33, p. 228], hence G is analytic on Sq. 
Obviously, we have <J>y ( (iy) = G(iy) for all y G M. d by Theorem 2.5 and the definition of G. Thus, [15, 
Lemma A.4] finally implies <J>y, = G on Sq- □ 
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Re 



o Js+ 



With Theorem 2.8 at hand, we can establish the following result: 

Theorem 2.9 (Absolute integrability). If (2.12) holds for some £ > 0, then w \-> <£>y, (y + iw) is abso- 
lutely integrable, for all y G M. d with \ \y\\ < B, where 6 is given as in Theorem 2.8. 

Proof. As in the proof of Theorem 2.8, we obtain from 

Re(H y+iw (s)) = H y (s) - I f 'e ATs ww T e As ds 

2 Jo 

and Re(> tr ( z )) < \e tT ^\ = e Rs ^ z » = e tr ^ z ^ for Z G M d (C), that 

M(H y+tw ( S ) + p*(y+i W ))X) _ A KL ( dX)d 

< J' J + ^W+P'MW-i) K L (dX)ds. 

Using this inequality yields 

\®Y t (y+iw)\ 

<^(y)e-2trto(«* T, A^(™^ 

= ^ r ( 3 ;)e-|(( A " 1 ^W(^o)+/oA- 1 ^W(u)^)w,H.) 

with SS[t) as in Lemma 2.4. Note that A~ l &(t)(Lo) + /„ A~ 1 ^(s)(y L )ds G Sj", hence 

\® Yt (y + iw)\dw < ® Y ,(y) [ e A{{^m^) + &x-^){ % )ds)^) dw < DOj 



by Theorem 2.8, and because the integrand is proportional to the density of a multivariate Normal 
distribution. □ 

2.5. Martingale Conditions and Equivalent Martingale Measures 

For notational convenience, we work in this section with the model 

dY, = (ji+Pfafidt+lfdWt+pidLt), Y Q £R d , (2.16) 
dL t = (Y L +AL t + L t A T )dt + dL t , L eS+ + , (2.17) 

where L is a driftless matrix subordinator with Levy measure Kl- Clearly, this is our multivariate 
stochastic volatility model of OU type (2.5), (2.6), except that p in (2.5) is replaced by p — p(/l), 
such that there is no deterministic drift from the leverage term p(dL t ). 

In mathematical finance, Y is used to model the joint dynamics of the log-returns of d assets with 
price processes S\ = S' Q e Y > , where we set Y Q ' = from now on and, hence, So denotes the vector of 
initial prices. 

The martingale property of the discounted stock prices (e~ rt S t ) te \oj] for a constant interest rate 
r > can be characterised as follows. 
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Theorem 2.10. The discounted price process [e rt S t \^)^\ is a martingale if and only if, for i = 
l,...,d, 

[ e p ' {x ^K L (dX) <oo, (2.18) 

J(\\x\\>i} 



and 



p i (X) = -^X ii , XGS+ (2.19) 



Ht = r- [ (eP'W-l)K L (dX). (2.20) 



Proof. Define 5, := e~ rt S t for all t £ [0, T] and let i G {1 , . . . , d}. By Ito's formula and [30, Proposition 
III.6.35], S' is a local martingale if and only if (2.18), (2.19) and (2.20) hold. Thus it remains to show 
that it is actually a true martingale under the stated assumptions. Since S is a positive local martingale, 
it is a supermartingale and hence a martingale if and only if E(S' T ) = Sq for all i G {l,...,d}. This 
can be seen as follows. By Theorem 2.3, (2.19) and (2.20) we have 

E(S' T ) = S' E (expfV-Or + ^E+H jfW)* dwj + p i W))) 



Sie~ Tk+ y (X) ~ X)KLm E ( e-mf+P'^E ( S 



se\0,r 



%e s d E[e p[LT >\ 



= S' Q 
= S' Q . 

This proves the assertion. □ 

As in [37, Theorem 3.1], it is possible to characterise the set of all equivalent martingale measures 
(henceforth EMMs), if the underlying filtration is generated by W and L. More specifically, it follows 
from the Martingale Representation Theorem (cf. [30, Theorem III.4.34]), that the density process 
Z t = E{^p\^ t ) of any equivalent martingale measure Q can be written as 

jT y s dW s + (F — 1) * Gu L - v L )^j (2.21) 

for suitable processes y and Y in this case. Here \i L resp. v L denote the random measure of jumps 
resp. its compensator (cf. [30, Section II. 1] for more details). Under an arbitrary EMM, L may not 
be a Levy process, and W and L may not be independent. However, there is a subclass of structure 
preserving EMMs under which L remains a Levy process independent of W. This translates into the 
following specifications of \j/ and Y (cf. [37, Theorem 3.2] for the univariate case): 

Theorem 2.11 (Structure preserving EMMs). Let y : St — > (0,°°) such that 



1- k +d {yW)-\) 2 KL{dX)<~., 

f(IWI>i} 



n 
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where k£(5) := f B y(X) K L (dX)forB£ 8(S$). Define the 
( 



Vt 



-E 2 



1 



V 



d -valued process (Yt)te[Q,T] as 

\ 



Ir 



J 



where 1 = (1, . . . , 1) T G W 1 . Then Z = S (J y s dW s + (y - 1) * (p L - V L )) is a density process, and 
the probability measure Q defined by ^ = Zj is an equivalent martingale measure. Moreover, 
W<2 := w — Jq \j/ s ds is a Q-standard Brownian motion, and L is an independent driftless Q-matrix 
subordinator with Levy measure tc[. The Q-dynamics of (Y, E) are given by 



dY l t 
dZ t 



1) Kl(dX) - Isf ) dt + ( z} dW t Q 



+ p i (dL t ), /=!,..., d, 



(YL+AL t + L t A T )dt + dL t . 



Proof. Since y— 1 > — l,Zis strictly positive by [30, Theorem 1.4.61]. The martingale property of Z 
follows along the lines of the proof of [37, Theorem 3.2]. The remaining assertions follow from [32, 
Proposition 1] and the Levy-Khintchine formula by applying the Girsanov-Jacod-Memin Theorem as 
in [32, Proposition 4] to the R2 d (^ +1 )-valued process 



W Q 




+ vech(L), 



where W Q := W - f y s ds. 



□ 



The previous theorem shows that it is possible to use a model of the same type under the real-world 
probability measure P and some EMM Q, e.g., to do option pricing and risk management within 
the same model class. The model parameters under Q can be determined by calibration, the model 
parameters under P by statistical methods. 



3. Option pricing using integral transform methods 

In this section we first recall results of [17] on Fourier pricing in general multivariate semimartingale 
models. To this end, let S = (5 ( )e y , . . . ,S^e Y ) be a J-dimensional semimartingale such that the dis- 
counted price process (e~ rt S t ) te [Qj] is a martingale under some pricing measure Q, for some constant 
instantaneous interest rate r > 0. 

We want to determine the price EQ(e~ rT f(Yj — s)) of a European option with payoff f{Yj — s) at 
maturity T, where / : W l — > M + is a measurable function and s := (— log(So), • • ■ > — l°§(^o))- Denote 
by / the Fourier transform of /. The following theorem is from [17, Theorem 3.2] and represents 
a multivariate generalisation of integral transform methods first introduced in the context of option 
pricing by [8] and [39]. 

Theorem 3.1 (Fourier Pricing). Fix R G M. d , let g(x) := e'^ f{x) for x G R d , and assume that 

(i) geL l D L°°, (ii) ®y t (R) < °°, (Hi) w H- <I>y r (R + iw) belongs to l> . 
Then, 

E Q (e- rT f(Y T -s)) = " ^ y J Rd e~'" M ®y t (R + iu)f(iR -u)du. (3.1) 
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Observe that Theorems 2.8 and 2.9 show that Conditions (ii) and (Hi) are satisfied for our multivari- 
ate stochastic volatility model of OU type (2.5), (2.6) if condition (2.12) holds, i.e., if L has enough 
exponential moments. More specifically, the vector R has to lie in the intersection of the domains of 
and /. 

We now present some examples. As is well-known, the Fourier transform of the payoff function of 
a plain vanilla call option with strike K > 0, f(x) = (e x — K) + is given by 

K x+iz 

f(z) = rj—r-, (3.2) 
k(1 + iz) 

for z £ C with Im(z) > 1. The Fourier transforms of many other single-asset options like barrier, 
self-quanto and power options as well as multi-asset options like worst-of and best-of options can be 
found, e.g., in the survey [17]. From the unpublished paper of [27] we have the following formulae 
for basket and spread options. 

Example 3.1. 1. The Fourier transform of fix) = (K — E/=i e *' ) + > K >0, that is the payoff func- 
tion of a basket put option, is given by 

u r(2+nf M zj) 

for all z £ C d with Im(zy) < 0, j = l,...,d. The price of the corresponding call can easily be 
derived using the put-call-parity (K — x) + = (x — K) + —x + K. Since we have separated the 
initial values s in (3.1), we can use F FT methods to compute the prices of weighted baskets for 
several weights efficiently. 

2. The Fourier transform of the payoff function of a spread call option, f(x) = (e* 1 — e Xl —K) + , 
K > 0, is given by 

= K l+iz,+iz2 T (iz 2 )T(-iz { -iZ2-l) 

nz) fei(i + fei) r(-fei-i) 

for all z £ C 2 with Im(zi) > 1, Im(z2) < and Im(zi +zi) > k see also [29]. 

Since the Fourier transform of (e X[ — e* 2 ) + does not exist anywhere, we cannot use Theorem 3.1 to 
price zero-strike spread options. Nevertheless, we can derive a similar formula directly. Alternatively, 
one could use the change of numeraire technique of [36], which would lead to formulae of a similar 
complexity. 

Proposition 3.2 (Spread options with zero strike). Suppose that 

*(y£,y£) {RA — R) < °° f or some R > 1 . 
Then the price of a zero-strike spread option with payoff (S^e^ — SqC Y t) + is given by 



„ , „ fP{si-si)-si-rT r .<b ( yi Y 2)(R + iuA—R — iu) 

E Q (e- rT (S T -S T ) + ) = - / g '»fe-.v,) r^-du, 

UK K T T! ' In Jr (R + iu)(R + iu-l) 

where s^ = — ln(5o) and S2 = — ln(5o). 

Observe that unlike for K > 0, one only has to compute a one-dimensional integral to determine the 
price of a zero-strike spread option. This will be exploited in the calibration procedure in Section 4. 
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Proof. Let R > 1 and define /k(x) = (e x — K) + for K > 0, and = e Rx /k(x). By Fourier inver- 
sion and (3.2), we have 

1 f e (R+iu)x e (l-R-iu)y 

f e y(x) = — -— — , w „ — : 77 du, 

w 2n.h (R + iu)(R + iu-\) 

for all }>GR. Hence, for the function h e y(x) := (S^e* — S$e y ) + = f e y-s 2 (x — si) we get 

h e y( x ) = -L^te-'i)-* / ,M&-*)-± -f -du. 

w 271 Jr (R + iu)(R + iu-\) 

Finally, by Fubini's theorem 

E Q (h ¥ 2 (Y})) = —— / / e^fe-'O * J _ 7sduP (Y i Y ^{dx,dy) 

u\ e y T \Tjj 2jl j r2 j r (R + iu)(R + iu-l) {YtJt > k ' 

'u) 

du, 



e R{s 2 -s l )-s 1 f <J>, yl y2 JR + iu, 1 — R - iu 

iu(s 2 —si) K't^t I v 



271 Jr (R + iu)(R + iu-l) 

where the application of Fubini's theorem is justified by 

f f e {R+iu)x e {\~R-iu)y 

/ / 7^ — • — 77T> — : 77 duP (Y i Y 2^dx,dy) = / e Rx e {l ~ R)y / \gi{u)\duP (Y i Y i,(dx,dy) 

J&h {R + iu)(R + iu-l) T t y> J R 2 V u n VtW » 

< \\gl\\v ®(Y},Yl)( R > l - R )<°°> 

since | \g\ | | L i < °° as shown in [17, Example 5. 1]. □ 

4. Calibration of the OU-Wishart model 

We now put forward a specific parametric specification of the model discussed in Section 2. To this 
end, let n £ N, £ §| and let X be a d x n random matrix with i.i.d. standard normal entries. Then, 
the matrix M := 05XX T 05 is said to be Wishart distributed, written M ~ Wd(n,@). Note that this 
definition can be extended to noninteger n > d — 1 using the characteristic function 

Zi-^det(^-2jZ0) _ 3" ) (4.1) 

see [23, Theorem 3.3.7]. Since M £ S J almost surely, we can define a compound Poisson matrix 
subordinator L with intensity X and Wd(n,&) distributed jumps. We call the resulting multivariate 
stochastic volatility model of OU type OU-Wishart model. 

Remark 4.1. There exists a subclass of structure preserving EMMs Q (cf. Theorem 2.11) such that 
we have an OU-Wishart model under both P and Q. This means that L is a compound Poisson 
process with Wd{n,©) distributed jumps and intensity X under P, and Wd(n,&) distributed jumps 
with intensity X under Q. We only need to assume that the Wishart distribution under both P and 
Q has a Lebesgue density, i.e., n,n > d — 1 and 0,0 £ Then, one simply has to take y as the 
quotient of the respective Levy densities. Hence, by [23, 3.2.1], y has to be defined as 

y(X) = I ( 2^ Fd Ufl det( ® } H ^ dei(X)^e-H^ 1 -^), Xg§|. 
A I M±n)det(0)W d 
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since we have / s+ e ir ^ K L {dX) = X det(/ d - 2R®)-<l" by (4.1), we see that the compound Poisson 

d 

process L has exponential moments as long as \\R\\ < jMjj' where ||-|| denotes the spectral norm. 

Consequently, (2.12) holds for e := jfpIT' anc ^ we can a PP^ tne integral transform methods from the 
previous section to compute prices of multi-asset options. 

Note that for the particularly simple special case of diagonal A, j3 and p, each asset follows a BNS 
model at the margins by (2.7) and (2.8). In particular, for n = 2 we see that L", i = l,...,d, is a 
compound Poisson subordinator with exponentially distributed jumps, thus we have in distribution 
the T-OU BNS model with stationary Gamma distribution at the margins, cf., e.g., [37, Section 2.2]. 
Then, the characteristic functions of the single assets are known in closed form. Note that while the 
characteristic function of the stationary distribution of the marginal OU type process is still known for 
n 7^ 2, it no longer corresponds to a Gamma distribution in this case. 



4.1. The OU-Wishart model in dimension 2 

We work directly under a pricing measure Q and consider the following specific two-dimensional case 
of our model, where we restrict ourselves in particular to a diagonal mean-reversion matrix A and a 
leverage term p such that both jumps of the respective variance and of the covariance enter the price. 
Our model is given by 

rfY/W/M l/^YU+fa 11 + (pidLV+pndL}* 

dY^-\W 2 )-2\n 2 )) dt+ W JfJ {dW 2 ) + [p 2 dI 22 + p 21 dLj 2 



dZ} 1 dZ] 2 \_((y l 0\ / 2a!!/ 1 ( fll +a 2 )Zj 2 \\ fdLj 1 dLj 2 
dZj 2 dl 22 ~ I I 7 2 J + W fll +a 2 )l/ 2 2a 2 Z 22 J J dt + \dL} 2 dl 22 



with initial values 



"0 

and parameters 71,72 > 0, a\,a 2 < 0, Pi,p 2 ,Pn,P2i 6i. L is a compound Poisson process with 
intensity X and ^2(«,0)-jumps, where n = 2 and 

V©12 ©22/ 

Therefore, all components of L jump at the same time. Since the second order properties of the 
Wishart distribution are known explicitly, cf. [23, Theorem 3.3.15], the covariances of the jumps are 
given by 

Cov(AL}\AL} 2 \AL? 7^0) =4011012, 
Cov(ALf ,AL/ 2 |AL/' ^0) =4022012, 
Cov ( AL/ 1 , ALj 2 1 AL/ 1 ^ 0) = 40 2 2 . 

This shows that even if p is diagonal, i.e., P12 = = P21, the leverage terms of both assets are con-el- 
ated. If p is non-diagonal, then di 2 also influences the marginal distribution of each asset. 
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Multi-asset option pricing By (2.14) and (4.1), the joint moment generating function of (Y l ,Y 2 ) 
is given by 



E[e> JY >] = 

exp Wilt + tr(Eofly(0) + f tr(nHy(s))ds + X jf 



1 



det(/ 2 -2(// v (s) + p*(y))0) 



ds — Xt 



with H y as in (2.15), A = ° J, y L = J), and p*(y) = (£* It does not seem to be 

possible to obtain a closed form expression in terms of ordinary functions, unless one sets a\ = a 2 =:a. 



In this case, if A = J 4b$b 2 — b\ ^ 0, one has 
E^V+yrf] = exp\y l p l t+y 2 p 2 t 4 



e 2at - 1 



4a 



tr Ei 



n-yi ym 
ym y\-y% 



+ T a -yi) + »0i -»)) (^(^' f - 1) ->) 



+ 



2abo 



— | arctan 
A 



2b 2 + b\ 



b + Z?i + Z? 2 
2 U^+V^ 



4m 



arctan 



+ — f — A? 



with coefficients 



&2 



l+4det(S-C) + 2tr(S-C), 

-8det(fl) +4tr(5)tr(C) - 4tr(BC) - 2tr(S), 

4det(S), 



and matrices 



5 



1 (y\-y\ ym 



4c? V ^1^2 yi -yi 



C:= (pm Pi 2 yA 
VP2U2 P2J2 / 



Note that arctan has to be understood as a function of complex argument to cover the case where the 
term in the square root of A is negative. If A = 0, we obtain 



exp \ y ip.it +y 2 il 2 t + 



? 2at _ l 



4a 



-tr L, 



n-yi ym 
ym y 2 -yi 



+T a (w(y?-yi) + 72 (A -»)) (^(^ - 1) -*) 



+ 



2a /3q 



ft 1 



2/3 2 e 2 "'+/ji 2b 2 + bi 2 \b 2 e 4at +b 1 e 2at + b 



bp + b\+b 2 

„4at 1 /,, ,>2a/ . 



+ — t-Xt 
bo 



Using det(A +B) = det(A) + det(5) +tr(A)tr(B) - tr(AB) for A,B E M 2 (R), the above formulae follow 
from 



det(/ 2 - 2(ffy(j) + p*(y))&) = det(/ 2 - 2(e 2ai -l)B-2C)=b + b x e las + /3 2 e to , 
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and straightforward integration. Likewise, one can also derive a closed form expression for n = 4, 6, . . . 
using [22,2.18(4)]. 

Consequently, one faces a tradeoff at this point. One possibility is to retain the flexibility of different 
mean reversion speeds a, by evaluating the remaining integral using numerical integration. Altern- 
atively, one can restrict attention to identical mean reversion speeds in order to have a closed-form 
expression of the moment generating function at hand. The impact of this decision on the calibration 
performance is discussed in Section 4.2 below. 

Single-asset option pricing For pricing single-asset options, one only needs the transforms of 
the marginal models, such that the above expressions simplify considerably. For example, the moment 
generating function of 7 1 is given by 

A , / b + b l \ Xt \ 

+ 2l^ ln l^T^J + ^- A 7' 

where bo and b\ simplify to 

^o = i+ ^2^Cyi-yi)-2piyi^ ©n-2pi 2 yi©i2, 
*i = -2^(y?-yi)eii. 

Note that one can use the recursion formula stated in [22, 2.155] to obtain a closed form expression 
for #2(«,@)-jumps with n € 2N, too. In the special case where the operator p is diagonal, i.e., 
if P12 = P21 = 0, the margins are (in distribution) T-OU BNS models, whose moment generating 
function has been derived in [37, Table 2.1]. 

Remark 4.2 (High Dimensionality). The above model can also be defined for d > 2, but of course, 
the Fourier formula (3.1) becomes numerically infeasible in high dimensions. Nevertheless, if p is 
diagonal, the calibration of a high dimensional OU-Wishart model is still possible by only evaluating 
options on two underlyings. Using zero strike spread options and provided the characteristic function 
is known explicitly, this means that one only has to evaluate single integrals numerically, as in the 
univariate case. Indeed, combining [4, Proposition 4.5] and the fact that every symmetric sub-matrix 
of a Wishart distributed matrix is again Wishart distributed, cf. [23, Theorem 3.3.10], it follows that 
the joint dynamics of each pair of assets follows a 2-dimensional OU-Wishart model as above. Hence, 
we can calibrate the model using only two-asset options (e.g., spread options). The price to pay is 
that the resulting model only incorporates pairwise dependencies, since the respective covariances 
completely determine the underlying Wishart distribution. 

Remark 4.3. If p is diagonal, we have equivalence in distribution of the margins of our model to a 
T-BNS model. This implies immediately that we need to use prices on multi-asset options in order to 
infer all parameters from observed option prices. If p is non-diagonal, we have a T-BNS model with 
an additional (correlated) jump term. Due to this additional term, it might be possible to infer ©12 
from single-asset options. However, one cannot obtain £q 2 in this way because it does not appear in 
the marginal moment generating function. 

In many multi-factor univariate models one can in general similarly not be sure whether one can 
uniquely determine all parameters from observed option prices. In many papers the parameters are 
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calibrated and the procedure seems to work, but we are not aware of any reasonably complex multi- 
factor model where the identifiability of the parameters based on option prices has been established. 
The reason is clearly the highly nontrivial relation between the parameters and the option prices. 

4.2. Empirical illustration 

The aim of this subsection is to show that a calibration of the OU-Wishart model to market prices 
is feasible. Since multi-asset options are mostly traded over-the-counter, it is difficult to obtain real 
price quotes. To circumvent this problem, we proceed as in [44] and consider foreign exchange rates 
instead, where a call option on some exchange rate can be seen as a spread option between two others. 
Let us emphasise that our calibration routine should not be seen as a finished product, but much 
rather as a first test and proof of principle. A more detailed investigation as well as an extension to 
numerically more involved models with non-diagonal A is left to future research. 

We consider a 2-dimensional OU-Wishart model as above. Our first asset is the EUR/USD ex- 
change rate S $//€ = S^/ € e Yl , that is, the price of 1 € in $, and our second asset is the GBP/USD 

exchange rate S $ / £ = S^J £ e Y ~ , i.e., the price of 1 £ in $. We model directly under a martingale meas- 
ure. Therefore we have, by Theorem 2. 10, that 

Ml = r $ - r € - [ (e* xU+ »* a - 1) K L (dX). 

Since Kl is the intensity X times a Wishart distribution with parameters n = 2 and 6, this simplifies to 

Ml =r$ - r € - X (det (l 2 - 2( P J ^ )0) 1 - 

2pi©n +2pi 2 012 

=r« — re — A, . 

* fe l-2p 1 0„-2p 12 12 

Likewise we have 

2p 2 2 2 + 2p 2 1012 



H2 = r$-r £ -X 



1 -2p 2 2 2 -2p 2 1012 



Thus, for P12 = or P21 = 0, we recover the martingale conditions of the T-OU BNS model. By 
[28, 13.4], it follows that the price in $ of a plain vanilla call option on S s / € or S s l £ is given by 
e- r $ T E((S $ T /€ - K)+) or e-'^ T E((S $ T /£ - K)+), respectively. Now observe that the $-payoff of a call 

option on the EUR/GBP exchange rate S £ / € is given by S $ T /£ (S £ T /€ - K)+ = (S $ T /€ - KS $ T /£ ) + , hence it 
can be regarded as a spread option on 5 $ / € — 5 $ / £ where the initial value of the second asset is replaced 

$/£ 

by KSfi . Since it is a zero-strike spread option, we can use Proposition 3.2 to valuate it. 

We obtained the option price data from EUWAX on April 29, 2010, at the end of the business 
day. The EUR/USD exchange rate at that time was sj§ /€ = 1.3249$, the GBP/USD exchange rate 

was sf/ £ = 1.5333$ and the EUR/GBP exchange rate was 0.8641£. As a proxy for the instantaneous 
riskless interest rate we took the 3-month LIBOR for each currency, viz. = 0.604%, r £ = 0.344% 
and r$ = 0.676%. All call options here are plain vanilla call options of European style. We used 148 
call options on the EUR/USD exchange rate, 67 call options on the GBP/USD exchange rate, and 105 
call options on the EUR/GBP exchange rate, all of them for different strikes and different maturities, 
for a total of 320 option prices. We always used the mid-value between bid and ask price. A spread 
sheet containing all data used for the calibration can be found on the second author's website. 
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ure 1: Comparison of the Black-Scholes implied volatility of market prices (dots) and model prices 
(solid line). The plots only show the results for the 12-parameter OU Wishart model (Step 
A), since they do not change visually for the more complex models from Step B to D. 
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Figure 2: Comparison of the Black-Scholes implied volatility of market prices (dots) and model prices 
(solid line). The plots only show the results for the 12-parameter OU Wishart model (Step 
A), since they do not change visually for the more complex models from Step B to D. 



The calibration was performed by choosing the model parameters so as to minimise the root mean 
squared error (RMSE) between the Black-Scholes volatilities implied by market resp. model prices. 
Note that the RMSE is the square root of the sum of the squared distances divided by the number of 
options. All computations were carried out in MATLAB and performed on a standard desktop PC 
with a 2.4GHz processor. 

In Step A, we impose a := a\ = aj and p\2 = = P21, i.e., we make the assumption that the mean 
reversion parameters of both assets are equal, and that p is diagonal. This is the most tractable case, 
since there is a closed form expression for the moment generating function of (Y 1 , Y 2 ) and the number 
of model parameters is reduced to 12. The starting and calibrated parameters can be found in Table 1. 
The overall RMSE is 0.0082, and the run time was 48 minutes, i.e., calibration of the model is feasible 
even on a standard PC. If one considers only the marginal models for EUR/USD and GBP/USD one 
has a RMSE of 0.0106 and 0.0048 respectively. For visualisation, we provide Figure 1 and 2, where 
market and model prices are compared in terms of Black-Scholes implied volatility for a few selected 
maturities. These results illustrate that even this simple model is able to fit the observed smiles rather 
well. For comparison, we calibrated two independent univariate T-OU BNS models to the margins 
separately (see Table 1) and obtained a lower RMSE of 0.0071 and 0.0020 respectively. This stems 
from the fact that the additional dependence parameters do not enter the pricing formulas for single 
asset options, whereas the intensity of the compound Poisson process is the same for all assets in 
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our multivariate framework, unlike when using two univariate models. This means that we are not 
overfitting the marginal distributions with an excessive amount of additional parameters, but much 
rather using a simplified version of a standard model. Nevertheless, the calibration still performs 
quite well even when using this simplification. 
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Figure 3: Simulated sample paths of the EUR/USD and the GBP/USD spot rates and their variances. 

As a further cross-check, Figure 3 depicts sample paths of the EUR/USD and the GBP/USD spot 
rates and their variances, simulated with our calibrated parameters, which show reasonable path prop- 
erties. 

In Step B, we allow for a non-diagonal leverage operator p . Although this introduces two additional 
parameters, p^ and P21 , a closed form expression for the moment generating function is still available. 
As initial values, we take the parameters obtained in Step A and set pu and P21 to zero. After 80 
minutes, the optimizer finds a minimum with a RMSE of 0.0079. At the margins, we have RMSEs 
of 0.0104 and 0.0037, respectively. Hence, calibration is still feasible without resorting to higher- 
powered computers, but the gains in fitting accuracy appear to be only moderate for the option price 
surface at hand. 

Next, we drop the assumption of an equal mean reversion parameter and allow for ci\ 7^ ai. Since 
the moment generating function of (Y l ,Y 2 ) is then not known in closed form anymore, good starting 
values are particularly important in order to reduce computational time to an acceptable value. We 
distinguish the two cases where p is diagonal (Step C) and p is non-diagonal (Step D), and take as 
starting values, the parameters obtained from Step A or Step B, respectively. Interestingly, in Step C 
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the optimizer finds the minimum at the same parameters as in Step A, thus the additional freedom of 
different mean reversion parameters does not yield a better fit in this case. 

Finally, in Step D, we calibrate the full model with non-diagonal p and different mean reversion 
speeds a\,a%. Due to the lack of a closed-form expression for the moment generating function and the 
high number of parameters (15), the run-time increases to an unsatisfactory 10 hours on our standard 
PC, suggesting that higher-powered computing facilities and an optimized numerical implementation 
in a compiled instead of an interpreted language should be employed here. In contrast to Step C, we 
find an improvement by allowing for different mean reversion speeds: The overall RMSE is 0.0076. 
Then again, for the data set at hand, the improvement is again only slight compared to the simplest 
model considered in Step A. 
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Table 1: Calibrated parameters for different models. In decreasing order: models from step A to D; 
univariate BNS model for EUR/USD and GBP/USD; initial parameters. 



Comparison With Other bivariate models We now compare our bivariate Wishart-OU model 
to some benchmarks from the literature. The canonical candidate would be the bivariate Wishart 
model, which also exhibits stochastic correlations between the assets and has very recently been cal- 
ibrated to market prices by [11]. However, the involved parameter restrictions necessary for the exist- 
ence of the Wishart process are not satisfied in the results of the calibration. This suggests that some 
kind of constrained optimization must be incorporated, which is beyond our scope here. However, we 
emphasize that the Wishart model should yield a comparable performance once these implementation 
issues have been resolved in a satisfactory manner. 

Instead, we use the multivariate Variance Gamma (henceforth VG) model of [34], and a generaliza- 
tion with stochastic volatility suggested therein for our comparison. In the mutivariate VG model with 
parameters (<?,-, <7,-, v), i = 1,2, the log -price processes Y l ,Y 2 are given by two independent Brownian 
motions with drift which are subordinated by a common Gamma process. The joint moment generat- 
ing function of the log-price processes under a risk neutral measure is shown to be given by 

E[^(y ^ Y t 1 +y2Y?)]=efrte- r ^ + y>te- r ^ , ( 1-vf (y^ + lyfaf 
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with wi = v _1 log (l — 6iV — |ofv). The parameters obtained from a calibration of this model to 
our option data set can be found in Table 2. The corresponding overall RMSE is 0.0134, which is 
roughly 63% higher than the RMSE obtained from the calibration of our 12-parameter OU-Wishart 
model from Step A. At the EUR/USD and GBP/USD margin the multivariate VG model has a RMSE 
of 0.0161 and 0.0107. Consequently, the performance of this model is much worse than for the OU- 
Wishart model, which is not surprising since it only involves 5 parameters. 

To alleviate this issue, our second benchmark allows for stochastic activity driven by an OU type 
process. More specifically, the log-price processes of the EUR/USD and GBP/USD spot rate are given 
by Y t = X\ t and Y t 2 = Xi, where X 1 and X are two independent Variance Gamma processes with 
parameters (0/, a,, v,-), i = 1,2, and Z t = J^z s ds is an integrated Ornstein-Uhlenbeck process. The 
Ornstein-Uhlenbeck process (z. v ) ie R+ is given by dz s = 2az s ds + dN-2at,Zo = 1, a < 0, where 7Y is a 
compound Poisson process with intensity ■& and Exp(<^) distributed jumps. It can be shown that the 
moment generating function of Z t , see, e.g., [42, 7.2.2], is given by 

= ex P ( f («p(2«) - 1) + 2g *fr - g/ °* [ - 2gg] iigfggg) - l Ji - 2 «w) . 

For the moment generating function of Y t = (Y t l ,Y t 2 ), conditioning on the stochastic activity process 
Z yields 

®y, (y i , yz ) = &z, (log ®xl (y i ) + log <& X 2 (y 2 ) ) 

with <E> Z j(y,) = (l — y,-0/Vj— \o}y\Vi) ^ V ', i = 1,2. Thus, the joint moment generating function of 
the log-price processes Y t l , Y t 2 under a risk neutral measure is given by 

^ Yt (i,oyy i ^ Yl (o,i)-^Y,(yi,y2). 
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Table 2: The first row shows the calibrated parameters for the multivariate VG model of [34]. The 
second row contains the calibrated parameters for two independent VG processes with a 
common integrated T-OU time change. 

A calibration of this model to our dataset leads to the parameters provided in Table 2; a plot depict- 
ing some of the respective implied volatilities can be found in Figure 4. The corresponding RMSE 
is 0.0129. Somewhat surprisingly, this is only around 4% lower than for the model of [34], despite 
increasing the parameters from 5 to 9. At the margins, we have 0.0143 and 0.0095, which corresponds 
to improvements of around 11%. Hence, there is quite some improvement in fitting the margins, but 
the multivariate options are not fit much better. This suggests that stochastic correlations indeed seem 
necessary to recapture the features of our empirical dataset. However, let us emphasize again that this 
only applies to one specific dataset in the foreign exchange market. A more detailed empirical study 
is a challenging topic for future research. 
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Figure 4: Comparison of the Black-Scholes implied volatility of market prices (dots) and model prices 
(solid line). The headers state the underlying and the days to maturity. The plots are for 
the benchmark model where the log-price processes are modelled by two independent VG 
processes with a common time change which is given by an integrated T-OU process. The 
plots for the multivariate VG model from [34] look very similar. 

5. Covariance swaps 

In this final section, we show that it is possible to price swaps on the covariance between different 
assets in closed form. This serves two purposes. On the one hand, options written on the realised co- 
variance represent a family of payoffs that only make sense in models where covariances are modeled 
as stochastic processes rather than constants. On the other hand, the ensuing calculations exemplify 
once more the analytical tractability of the present framework. 

We consider again our multivariate stochastic volatility model of OU type under an EMM Q. In ad- 
dition, we suppose that the matrix subordinator L is square integrable, i.e., f^ x \\>\} \\^-\\ 2 K L(dX) < oo. 
The pricing of options written on the realised variance resp. the quadratic variation as its continuous- 
time limit have been studied extensively in the literature, cf., e.g., [6] and the references therein. 
Since we have a nontrivial correlation structure in our model, one can also consider covariance swaps 
on two assets i,j G {l,...,d}, i.e., contracts with payoff [F^F^y — K with covariance swap rate 
K = E{\Y\Y^\t) (see, e.g., [7], [12], or [43] for more background on these products). Now, we show 
how to compute the covariance swap rate. We have 

[F',W]r = rjJfr + £ AY'AYj = , + p i {X)p\X) *${dX). 

s<T 
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Since Ki{dX)dt is the compensator of \i L , this yields 

= (£(£+)) ! 'J + T I p i (X)pJ(X)K L (dX), 



(5.1) 



where Lj was defined in Equation (2.4). Note that by [38, Proposition 2.4] and since \p'(X)p->(X)\ < 
| \p\ \ 2 \\X\ | 2 , our integrability assumption on L implies that the expectation is finite. The first summand 
can be calculated as follows. By setting y = in Theorem 2.5, we obtain the characteristic function 
of E f . Differentiation yields 

E(-L T ) = <FlZ Q <? r + e AT A-\E(L l ))e AJT - A" 1 (E(Li)), 

where E{L\) = Jl + J§+X Ki(dX). Using Equation (2.4), we obtain 

£(£+)= A-\E(L T ) - TE{L\) - Eo), 

so we only need to know E{L\). The second summand in (5.1) can analogously be computed by 
differentiating the characteristic function of the matrix subordinator L. 

In our OU-Wishart model, where L is a compound Poisson matrix subordinator plus drift with 
Wd(n,®) -distributed jumps, we have by [23, Theorem 3.3.15] that 

E(L l ) = y L + Xn&. 

If p is diagonal, the second term in (5.1) simplifies to 

Tpipj f ^XnXjjVidX) = TpiPjXn (20?. + n© ii © jj ) , 
again by [23, Theorem 3.3.15]. Thus we have a closed form expression for the covariance swap rate: 

K = (A- 1 \e AT {L Q + A-\y L + Xn<d))e ATT -A-\y L + Xn<d)-T{y L + Xn<d)-^ ,} 
+ Tp i pjXn(2&l + n& ii & jj ) . 

For example, in the 2-dimensional OU-Wishart model from Section 4.1 we have, for i = 1 and j = 2, 
1 



K 



U*i+°i)T_A (^ + ^®n\_ TXn&i2 \ + Tp l p 2 Xn(2® 2 l2 + n® n ®22) 



CL\ +C12 

As an illustration we provide, in Figure 5, a plot of the normalized covariance swap rate measured in 



volaility points, i.e., T i->- y jE([Y 1 ,Y 2 ]t), for our calibrated 12-parameter OU-Wishart model from 
Section 4.2 (Step A). 

Finally, we remark that similarly as in [6], pricing of options on the covariance can be dealt with 
using the Fourier methods from Section 3, since the joint characteristic function of (L + ,p'(X)pj(X) * 
p L {dX)) can be calculated similarly as in the proof of Theorem 2.5. 
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Normalized Covariance Swap Rate in Voaltility Points 
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Figure 5: Normalized covariance swap rate for the calibrated 12-parameter OU-Wishart model. 

A. Appendix 

The following result on multidimensional analytic functions is needed in the proof of Lemma 2.7. 

Lemma A.l. Let D £ = {z G C" : | |Re(z)|| < e} for some e > 0. Suppose f : D e — >■ C is an analytic 
function of the form f = e F , where F : D e — > C is continuous. Then F is analytic in D e . 

Proof. Let z = {zi,Zi, ■ ■ ■ ,z n ) € T) e and define z_i = (z2,- ■ ■ ,z„). Then / z _, : w i-> /(w,z_i) defines 
an analytic function without zeros on the open convex set D £ z _, := {w 6 C : (w,z_i) € O e }. By, e.g., 
[19, Satz V.1.4], there exists an analytic function : D E z { — > C such that e\p(gl_ ) = f z _ r Hence 
F(w,z~i) ~8z-i ( w ) e 27T/Z on D e z r Since both F and g are continuous, their difference is constant 
and it follows that w i-)- F(w,z_i) is analytic on D e z _ r Analogously, one shows analyticity of F in 
all other components. The assertion then follows from Hartog's Theorem (cf., e.g., [25, Theorem 
2.2.8]). □ 
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